Tridiagonal matrix algorithm

 ( Numerical Linear Algebra ) 

$c\text{'}\_i\; =$

\begin{cases}
 \cfrac{c_i}{b_i},                  & i = 1, \\
 \cfrac{c_i}{b_i - a_i c'_{i - 1}}, & i = 2, 3, \dots, n - 1
\end{cases}
     

and

$d\text{'}\_i\; =$

\begin{cases}
 \cfrac{d_i}{b_i},                                   & i = 1, \\
 \cfrac{d_i - a_i d'_{i - 1}}{b_i - a_i c'_{i - 1}}, & i = 2, 3, \dots, n.
\end{cases}
     

The solution is then obtained by back substitution:

$x\_n\; =\; d\text{'}\_n,$

$x\_i\; =\; d\text{'}\_i\; -\; c\text{'}\_i\; x\_\{i\; +\; 1\},\; \backslash quad\; i\; =\; n\; -\; 1,\; n\; -\; 2,\; \backslash ldots,\; 1.$

The method above does not modify the original coefficient vectors, but must also keep track of the new coefficients. If the coefficient vectors may be modified, then an algorithm with less bookkeeping is:

For $i\; =\; 2,\; 3,\; \backslash dots,\; n,$ do

$w\; =\; \backslash cfrac\{a\_i\}\{b\_\{i-1\}\},$

$b\_i\; :=\; b\_i-w\; c\_\{i-1\},$

$d\_i\; :=\; d\_i-w\; d\_\{i-1\},$

followed by the back substitution

$x\_n\; =\; \backslash cfrac\{d\_n\}\{b\_n\},$

$x\_i\; =\; \backslash cfrac\{d\_i\; -\; c\_i\; x\_\{i+1\}\}\{b\_i\}\; \backslash quad\; \backslash text\{for\; \}\; i\; =\; n\; -\; 1,\; n\; -\; 2,\; \backslash dots,\; 1.$

The implementation as a C function, which uses scratch space to avoid modifying its inputs for a-c, allowing them to be reused: void thomas(const int X, double xrestrict,

           const double a[restrict X], const double b[restrict X],
           const double c[restrict X], double scratch[restrict X]) {
   /*
    solves Ax = d, where A is a tridiagonal matrix consisting of vectors a, b, c
    X = number of equations
    x[] = initially contains the input, d, and returns x. indexed from [0, ..., X - 1]
    a[] = subdiagonal, indexed from [1, ..., X - 1]
    b[] = main diagonal, indexed from [0, ..., X - 1]
    c[] = superdiagonal, indexed from [0, ..., X - 2]
    scratch[] = scratch space of length X, provided by caller, allowing a, b, c to be const
    not performed in this example: manual expensive common subexpression elimination
    */
   scratch[0] = c[0] / b[0];
   x[0] = x[0] / b[0];
     

   /* loop from 1 to X - 1 inclusive */
   for (int ix = 1; ix < X; ix++) {
       if (ix < X-1){
       scratch[ix] = c[ix] / (b[ix] - a[ix] * scratch[ix - 1]);
       }
       x[ix] = (x[ix] - a[ix] * x[ix - 1]) / (b[ix] - a[ix] * scratch[ix - 1]);
   }
     

   /* loop from X - 2 to 0 inclusive */
   for (int ix = X - 2; ix >= 0; ix--)
       x[ix] -= scratch[ix] * x[ix + 1];
     

Suppose that the unknowns are $x\_1,\backslash ldots,\; x\_n$, and that the equations to be solved are:

$\backslash begin\{alignat\}\{4\}$

Consider modifying the second ($i\; =\; 2$) equation with the first equation as follows:

which would give:

Note that $x\_1$ has been eliminated from the second equation. Using a similar tactic with the modified second equation on the third equation yields:

= d_3 (b_2 b_1 - c_1 a_2) - (d_2 b_1 - d_1 a_2) a_3.
     

This time $x\_2$ was eliminated. If this procedure is repeated until the $n^\{th\}$ row; the (modified) $n^\{th\}$ equation will involve only one unknown, $x\_n$. This may be solved for and then used to solve the $(n\; -\; 1)^\{th\}$ equation, and so on until all of the unknowns are solved for.

Clearly, the coefficients on the modified equations get more and more complicated if stated explicitly. By examining the procedure, the modified coefficients (notated with tildes) may instead be defined recursively:

$\backslash tilde\; a\_i\; =\; 0\backslash ,$

$\backslash tilde\; b\_1\; =\; b\_1\backslash ,$

$\backslash tilde\; b\_i\; =\; b\_i\; \backslash tilde\; b\_\{i\; -\; 1\}\; -\; \backslash tilde\; c\_\{i\; -\; 1\}\; a\_i\backslash ,$

$\backslash tilde\; c\_1\; =\; c\_1\backslash ,$

$\backslash tilde\; c\_i\; =\; c\_i\; \backslash tilde\; b\_\{i\; -\; 1\}\backslash ,$

$\backslash tilde\; d\_1\; =\; d\_1\backslash ,$

$\backslash tilde\; d\_i\; =\; d\_i\; \backslash tilde\; b\_\{i\; -\; 1\}\; -\; \backslash tilde\; d\_\{i\; -\; 1\}\; a\_i.\backslash ,$

To further hasten the solution process, $\backslash tilde\; b\_i$ may be divided out (if there's no division by zero risk), the newer modified coefficients, each notated with a prime, will be:

$a\text{'}\_i\; =\; 0\backslash ,$

$b\text{'}\_i\; =\; 1\backslash ,$

$c\text{'}\_1\; =\; \backslash frac\{c\_1\}\{b\_1\}\backslash ,$

$c\text{'}\_i\; =\; \backslash frac\{c\_i\}\{b\_i\; -\; c\text{'}\_\{i\; -\; 1\}\; a\_i\}\backslash ,$

$d\text{'}\_1\; =\; \backslash frac\{d\_1\}\{b\_1\}\backslash ,$

$d\text{'}\_i\; =\; \backslash frac\{d\_i\; -\; d\text{'}\_\{i\; -\; 1\}\; a\_i\}\{b\_i\; -\; c\text{'}\_\{i\; -\; 1\}\; a\_i\}.\backslash ,$

This gives the following system with the same unknowns and coefficients defined in terms of the original ones above:

$\backslash begin\{array\}\{lcl\}$

The last equation involves only one unknown. Solving it in turn reduces the next last equation to one unknown, so that this backward substitution can be used to find all of the unknowns:

$x\_n\; =\; d\text{'}\_n\backslash ,$

$x\_i\; =\; d\text{'}\_i\; -\; c\text{'}\_i\; x\_\{i\; +\; 1\}\; \backslash qquad\; ;\; \backslash \; i\; =\; n\; -\; 1,\; n\; -\; 2,\; \backslash ldots,\; 1.$

In this case, we can make use of the Sherman–Morrison formula to avoid the additional operations of Gaussian elimination and still use the Thomas algorithm. The method requires solving a modified non-cyclic version of the system for both the input and a sparse corrective vector, and then combining the solutions. This can be done efficiently if both solutions are computed at once, as the forward portion of the pure tridiagonal matrix algorithm can be shared.

If we indicate by: $$A=\backslash begin\{bmatrix\}$$

  b_1 & c_1 &        &        &  a_1      \\
  a_2 & b_2 & c_2    &        &         \\
      & a_3 & b_3    & \ddots &         \\
      &     & \ddots & \ddots & c_{n-1} \\
  c_n   &     &        & a_n    & b_n
     

  x_1    \\
  x_2    \\
  x_3    \\
  \vdots \\
  x_n
     

  d_1    \\
  d_2    \\
  d_3    \\
  \vdots \\
  d_n
     

Then the system to be solved is: $$Ax\; =\; d$$

In this case the coefficients $a\_1$ and $c\_n$ are, generally speaking, non-zero, so their presence does not allow to apply the Thomas algorithm directly. We can therefore consider $B\backslash in\backslash mathbb\{R\}^\{n\backslash times\; n\}$ and $u,v\backslash in\backslash mathbb\{R\}^\{n\}$ as following: $$B=\backslash begin\{bmatrix\}$$

  b_1-\gamma & c_1 &        &        &  0      \\
  a_2 & b_2 & c_2    &        &         \\
      & a_3 & b_3    & \ddots &         \\
      &     & \ddots & \ddots & c_{n-1} \\
  0   &     &        & a_n    & b_n-\frac{c_na_1}{\gamma}
     

  \gamma    \\
  0    \\
  0    \\
  \vdots \\
  c_n
     

  1    \\
  0    \\
  0    \\
  \vdots \\
  a_1/\gamma
     

Then we reconstruct the solution using the Shermann-Morrison formula:

The implementation as a C function, which uses scratch space to avoid modifying its inputs for a-c, allowing them to be reused: void cyclic_thomas(const int X, double xrestrict, const double arestrict, const double brestrict, const double crestrict, double cmodrestrict, double urestrict) {

   /*
    solves Ax = v, where A is a cyclic tridiagonal matrix consisting of vectors a, b, c
    X = number of equations
    x[] = initially contains the input v, and returns x. indexed from [0, ..., X - 1]
    a[] = subdiagonal, regularly indexed from [1, ..., X - 1], a[0] is lower left corner
    b[] = main diagonal, indexed from [0, ..., X - 1]
    c[] = superdiagonal, regularly indexed from [0, ..., X - 2], c[X - 1] is upper right
    cmod[], u[] = scratch vectors each of length X
    */
     

   /* lower left and upper right corners of the cyclic tridiagonal system respectively */
   const double alpha = a[0];
   const double beta = c[X - 1];
     

   /* arbitrary, but chosen such that division by zero is avoided */
   const double gamma = -b[0];
     

   cmod[0] = c[0] / (b[0] - gamma);
   u[0] = gamma / (b[0] - gamma);
   x[0] /= (b[0] - gamma);
     

   /* loop from 1 to X - 2 inclusive */
   for (int ix = 1; ix + 1 < X; ix++) {
       const double m = 1.0 / (b[ix] - a[ix] * cmod[ix - 1]);
       cmod[ix] = c[ix] * m;
       u[ix] = (0.0f  - a[ix] * u[ix - 1]) * m;
       x[ix] = (x[ix] - a[ix] * x[ix - 1]) * m;
   }
     

   /* handle X - 1 */
   const double m = 1.0 / (b[X - 1] - alpha * beta / gamma - a[X - 1] * cmod[X - 2]);
   u[X - 1] = (alpha    - a[X - 1] * u[X - 2]) * m;
   x[X - 1] = (x[X - 1] - a[X - 1] * x[X - 2]) * m;
     

   /* loop from X - 2 to 0 inclusive */
   for (int ix = X - 2; ix >= 0; ix--) {
       u[ix] -= cmod[ix] * u[ix + 1];
       x[ix] -= cmod[ix] * x[ix + 1];
   }
     

   const double fact = (x[0] + x[X - 1] * alpha / gamma) / (1.0 + u[0] + u[X - 1] * alpha / gamma);
     

   /* loop from 0 to X - 1 inclusive */
   for (int ix = 0; ix < X; ix++)
       x[ix] -= fact * u[ix];
     

There is also another way to solve the slightly perturbed form of the tridiagonal system considered above. Let us consider two auxiliary linear systems of dimension $(n-1)\; \backslash times\; (n-1)$: $$\backslash begin\{align\}$$

\qquad \ \ \ \ \ b_2 u_{2}  + c_2 u_3 &= d_2 \\
a_3 u_2 + b_3 u_3 + c_3 u_4  &= d_3 \\
a_i u_{i-1} + b_i u_i + c_i u_{i+1} &= d_i\\
\dots \\
a_n u_{n - 1}+ b_n u_n\qquad &= d_n \,.
     

\begin{align}

\qquad \ \ \ \ \ b_2 v_{2}  + c_2 v_3 &= -a_2 \\
a_3 v_2 + b_3 v_3 + c_3 v_4  &= 0 \\
a_i u_{i-1} + b_i u_i + c_i u_{i+1} &= 0\\
\dots \\
a_n v_{n - 1}+ b_n v_n\qquad &= -c_n \,.
     

For convenience, we additionally define $u\_1\; =\; 0$ and $v\_1\; =\; 1$. We can now find the solutions $\backslash \{u\_2,u\_3\backslash dots,u\_n\; \backslash \}$ and $\backslash \{v\_2,v\_3\backslash dots,v\_n\; \backslash \}$ applying Thomas algorithm to the two auxiliary tridiagonal system.

The solution $\backslash \{\; x\_1,x\_2\backslash dots,x\_n\; \backslash \}$ can be then represented in the form: $$x\_i\; =\; u\_i\; +\; x\_1\; v\_i\; \backslash qquad\; i=1,2,\backslash dots,\; n$$

Indeed, multiplying each equation of the second auxiliary system by $x\_1$, adding with the corresponding equation of the first auxiliary system and using the representation $x\_i\; =\; u\_i\; +\; x\_1\; v\_i$, we immediately see that equations number through of the original system are satisfied; it only remains to satisfy equation number . To do so, consider formula for $i=2$ and $i=n$ and substitute $x\_2\; =\; u\_2\; +\; x\_1\; v\_2$and $x\_n\; =\; u\_n\; +\; x\_1\; v\_n$ into the first equation of the original system. This yields one scalar equation for $x\_1$: $$b\_1x\_1+c\_1(u\_2+x\_1v\_2)+a\_1(u\_n+x\_1v\_n)\; =\; d\_1$$

As such, we find: $$x\_1\; =\; \backslash frac\{d\_1-a\_1u\_n-c\_1u\_2\}\{b\_1+a\_1v\_n+c\_1v\_2\}$$

The implementation as a C function, which uses scratch space to avoid modifying its inputs for a-c, allowing them to be reused: void cyclic_thomas(const int X, double xrestrict, const double arestrict, const double brestrict, const double crestrict, double cmodrestrict, double vrestrict) {

   /* first solve a system of length X - 1 for two right hand sides, ignoring ix == 0 */
   cmod[1] = c[1] / b[1];
   v[1] = -a[1] / b[1];
   x[1] = x[1] / b[1];
     

   /* loop from 2 to X - 1 inclusive */
   for (int ix = 2; ix < X - 1; ix++) {
       const double m = 1.0 / (b[ix] - a[ix] * cmod[ix - 1]);
       cmod[ix] = c[ix] * m;
       v[ix] = (0.0f  - a[ix] * v[ix - 1]) * m;
       x[ix] = (x[ix] - a[ix] * x[ix - 1]) * m;
   }
     

   /* handle X - 1 */
   const double m = 1.0 / (b[X - 1] - a[X - 1] * cmod[X - 2]);
   cmod[X - 1] = c[X - 1] * m;
   v[X - 1] = (-c[0]    - a[X - 1] * v[X - 2]) * m;
   x[X - 1] = (x[X - 1] - a[X - 1] * x[X - 2]) * m;
     

   /* loop from X - 2 to 1 inclusive */
   for (int ix = X - 2; ix >= 1; ix--) {
       v[ix] -= cmod[ix] * v[ix + 1];
       x[ix] -= cmod[ix] * x[ix + 1];
   }
     

   x[0] = (x[0] - a[0] * x[X - 1] - c[0] * x[1]) / (b[0] + a[0] * v[X - 1] + c[0] * v[1]);
     

   /* loop from 1 to X - 1 inclusive */
   for (int ix = 1; ix < X; ix++)
       x[ix] += x[0] * v[ix];
     

In both cases the auxiliary systems to be solved are genuinely tri-diagonal, so the overall computational complexity of solving system $Ax\; =\; d$ remains linear with the respect to the dimension of the system , that is $O(n)$ arithmetic operations.

In other situations, the system of equations may be block tridiagonal (see block matrix), with smaller submatrices arranged as the individual elements in the above matrix system (e.g., the 2D Poisson problem). Simplified forms of Gaussian elimination have been developed for these situations.